All sky research

Solo ∼ 3 × 10 3 stelle di neutroni osservate su ∼ 10 9 stimate nella Via Lattea ⇒ ricerca “alla cieca” su tutto il cielo. Problemi

  • ampiezze molto piccole;
  • modulazione Doppler dovuto ai moti terrestri;
  • dimensione dello spazio dei parametri: frequenza, spin-down, posizione nel cielo.

Frequency-Hough pipeline

  • Trasformate di Fourier discrete tali che ∆ν < ∆ν$_{Doppler}$ (8192 s, 4096 s)
  • Mappa tempo-frequenza (peakmap)
  • Correzione Doppler per ogni posizione del cielo
  • Trasformata Frequency-Hough
  • Selezione candidati

Astone, P. et al., Phys. Rev. D 90.4, 042002 (2014).

Il codice originale


In [ ]:
binh_df0ORIG=zeros(nbin_d,nbin_f0);  %  HM matrix container
for it = 1:nTimeSteps
    kf=(peaks(2,ii0:ii(it))-inifr)/ddf;  % normalized frequencies
    w=peaks(5,ii0:ii(it));               % wiener weights
    t=peaks(1,ii0)*Day_inSeconds; % time conversion days to s
    tddf=t/ddf;
    f0_a=kf-deltaf2; 
    for id = 1:nbin_d   % loop for the creation of half-differential map
        td=d(id)*tddf;
        a=1+round(f0_a-td+I500); 
        binh_df0ORIG(id,a)=binh_df0ORIG(id,a)+w; % left edge of the strips
    end
    ii0=ii(it)+1;
end

binh_df0ORIG(:,deltaf2*2+1:nbin_f0)=...
    binh_df0ORIG(:,deltaf2*2+1:nbin_f0)-binh_df0ORIG(:,1:nbin_f0-deltaf2*2); % half to full diff. map - Carl Sabottke idea
binh_df0ORIG=cumsum(binh_df0ORIG,2);   % creation of the Hough map

Il nuovo codice


In [ ]:
def rowTransform(ithSD):
    sdTimed = tf.multiply(spindowns[ithSD], times)
    transform = tf.round(frequencies-sdTimed+securbelt/2)
    transform = tf.cast(transform, dtype=tf.int32)
    values = tf.unsorted_segment_sum(weights, transform, numColumns)
    values = tf.cast(values, dtype=tf.float32)
    return values

houghLeft = tf.map_fn(rowTransform, tf.range(0, numRows), dtype=tf.float32, parallel_iterations=8)
houghRight = houghLeft[:,enhancement:numColumns]-houghLeft[:,0:numColumns - enhancement]
houghDiff = tf.concat([houghLeft[:,0:enhancement],houghRight],1)
houghMap = tf.cumsum(houghDiff, axis = 1)

Let's review it


In [ ]:
#this function computes the Hough transform histogram for a given spindown
def rowTransform(ithSD):
    sdTimed = tf.multiply(spindowns[ithSD], times)

    transform = tf.round(frequencies-sdTimed+securbelt/2)
    transform = tf.cast(transform, dtype=tf.int32)
    #the rounding operation brings a certain number of peaks in the same 
    #frequency-spindown bin in the Hough map
    #the left edge is then computed binning that peaks properly
    #(according to their values if the peakmap was adactive)
    #this is the core of the algoritm and brings the most computational effort
    values = tf.unsorted_segment_sum(weights, transform, numColumns)
    values = tf.cast(values, dtype=tf.float32)
    return values
#to keep under control the memory usage, the map function is a 
#very useful tool to apply the same function over a vector
#in this way the vectorization is preserved
houghLeft = tf.map_fn(rowTransform, tf.range(0, numRows), 
                      dtype=tf.float32, parallel_iterations=8)
#let's superimpose the right edge on the image
houghRight = houghLeft[:,enhancement:numColumns]-houghLeft[:,0:numColumns - enhancement]
houghDiff = tf.concat([houghLeft[:,0:enhancement],houghRight],1)
#at last, the Hough map is computed integrating along the frequencies
houghMap = tf.cumsum(houghDiff, axis = 1)

How it works